Reproducible research involves analysis and reporting of the data from raw data to final outcomes in one unified framework. Given the data, another researcher should be able to run it through the same analysis, and see each step of the analysis and check the robustness of the findings as they wish. This also means that the addition of extra data, or omission of problematic data, or changes to the way that data should be processed should be straightforward to implement without the need of error-prone manual entering and reporting of data. In this way, research is made more transparent, errors should be rarer and easier to diagnose, and the analysis itself is also made easier in the long run, with its requiring little to no extra work as data sets increase, or errors in the raw data are discovered.
R is the ideal option for performing reproducible research as its primary IDE, RStudio, comes prepackaged with R Markdown and knitr, which is specifically designed for reproducible research and generating reproducible reports. Furthermore, as opposed to MATLAB, it is free which means that the data analysis can be run by anyone and not only those who pay for the software. It is also a language designed for statistical analysis, meaning that it is ideally suited to the statistical analysis and presentation. Finally, there are powerful tools within R for performing the fitting procedures required for effective kinetic modelling.
Software available at present for TAC extraction and kinetic analysis include commercial products, MIAKAT (which is open-source, but still requires MATLAB which is not freeware), and various custom-built tools within research groups, most often built in MATLAB. Statistical modelling and presentation are sometimes performed in MATLAB, but are usually performed using various statistical and graphical programs.
All the processes up until TAC extraction require large raw data files to be shared and to be reproducible, which are often complicated not only to store, but also to share as they contain sensitive patient information in the form of an image, which is theoretically difficult to fully anonymise and therefore ethically problematic. TAC extraction is also a process which can be performed in many different software packages as this does not differ between PET and fMRI, including, for example, NiPy in Python, which is also a free software package. Storage and sharing of TACs is more convenient as this data is in the form of a text document containing one-dimensional time series data. This therefore, in the opinion of the author, represents the ideal step of the image analysis pipeline to begin a reproducible analysis which can be shared in full. This can also be kept for documentation purposes in the interests of data transparency.
For these reasons, and to enhance my own knowledge of popular kinetic modelling procedures, I have written a kinetic modelling R package which can analyse and plot TAC data using several popular kinetic models, while easily accepting input data for each function in as raw a form as possible without requiring complex data structures. Furthermore, each process stores not only the parameters, but also the fit itself, the fitted values, all of the input data, and has straightforward tools for plotting the model to perform quality control.
# Loading packages
library(kinfitr)
library(knitr)
library(pander)
library(broom)times = tacdata$Time/60
reftac <- tacdata$RefCBL
roitac <- tacdata$FSLSSTR
weights <- tacdata$WeightsFirst all the data is loaded (unshown). This data includes:
This data is plotted below:
par(mfrow=c(3,1))
plot(times, reftac)
plot(times, roitac)
plot(times, weights)Reference Region Model Data
\[C(t) = R_{1}C_{R}(t)+(k_{2}-\frac{R_{1}k_{2}}{1+BP_{ND}})C_{R}(t)\otimes e^{-\frac{k_{2}}{1+BP_{ND}}t}\]
SRTM requires input of the times, reftac and roitac. Weights are optional.
srtmout <- srtm(times, reftac, roitac, weights)pandoc.table(srtmout$par, digits = 3)| R1 | k2 | bp |
|---|---|---|
| 0.856 | 0.111 | 1.53 |
The output object can be used to visualise the fit:
plot_srtmfit(srtmout)SRTM Fit Check
… and the residuals:
plot_residuals(srtmout)SRTM Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(srtmout))| Length | Class | Mode | |
|---|---|---|---|
| par | 3 | data.frame | list |
| par.se | 3 | data.frame | list |
| fit | 6 | nls | list |
| weights | 14 | -none- | numeric |
| tacs | 4 | data.frame | list |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(srtmout$fit), digits=3)| sigma | isConv | finTol | logLik | AIC | BIC | deviance | df.residual |
|---|---|---|---|---|---|---|---|
| 5.43 | TRUE | 1.49e-08 | -57.4 | 123 | 125 | 295 | 10 |
refLogan requires input of the times, reftac, roitac and k2prime. Weights are optional.
k2prime can be obtained from the srtm fit:
k2prime = srtmout$par$k2 / srtmout$par$R1
k2prime## [1] 0.1292453
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
taclow = tacdata$FSLSFC # Frontal Cortex
tacmed = tacdata$fslLST # Limbic Striatum
tachigh = tacdata$fslAST # Associative Striatum
refLogan_tstar(times, reftac, taclow, tacmed, tachigh, k2prime = k2prime, 'demonstration')## png
## 2
Using this information, we see that we should probably use 9 included frames for the fit.
refloganout <- refLogan(times, reftac, roitac, k2prime,
tstarIncludedFrames = 9, weights=weights)pandoc.table(refloganout$par)| bp |
|---|
| 1.522 |
The output object can be used to visualise the fit:
plot_refLoganfit(refloganout)RefLogan Fit Check
… and the residuals:
plot_residuals(refloganout)RefLogan Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(refloganout))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 3 | data.frame | list |
| fitvals | 2 | data.frame | list |
| weights | 14 | -none- | numeric |
| k2prime | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(refloganout$fit), digits = 3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 0.279 | 40528 | 1.97e-14 | 2 | -0.742 | 7.48 |
| BIC | deviance | df.residual |
|---|---|---|
| 8.08 | 0.543 | 7 |
MRTM1 requires input of the times, reftac and roitac. Weights are optional.
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
mrtm1_tstar(times, reftac, taclow, tacmed, tachigh, 'demonstration')Using this information (unshown), we see that we should probably use all 14 frames for the fit.
mrtm1out <- mrtm1(times, reftac, roitac, tstarIncludedFrames = 14, weights=weights)pandoc.table(mrtm1out$par, digits = 3)| bp | k2prime |
|---|---|
| 1.52 | 0.13 |
The output object can be used to visualise the fit:
plot_mrtm1fit(mrtm1out)MRTM1 Fit Check
… and the residuals:
plot_residuals(mrtm1out)MRTM1 Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(mrtm1out))| Length | Class | Mode | |
|---|---|---|---|
| par | 2 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 3 | data.frame | list |
| fitvals | 9 | data.frame | list |
| weights | 14 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(mrtm1out$fit), digits = 3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 5.1 | 16740 | 8.46e-19 | 3 | -39.2 | 86.4 |
| BIC | deviance | df.residual |
|---|---|---|
| 88.7 | 260 | 10 |
MRTM2 requires input of the times, reftac, roitac and k2prime. Weights are optional.
k2prime can be obtained from the mrtm1 fit:
k2prime = mrtm1out$par$k2prime
k2prime## [1] 0.1295072
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
mrtm2_tstar(times, reftac, taclow, tacmed, tachigh, k2prime=k2prime, 'demonstration')Using this information (unshown), we see that we should probably use all 14 frames for the fit.
mrtm2out <- mrtm2(times, reftac, roitac, k2prime = k2prime,
tstarIncludedFrames = 14, weights=weights)pandoc.table(mrtm2out$par)| bp |
|---|
| 1.521 |
The output object can be used to visualise the fit:
plot_mrtm2fit(mrtm2out)MRTM2 Fit Check
… and the residuals:
plot_residuals(mrtm2out)MRTM2 Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(mrtm2out))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 3 | data.frame | list |
| fitvals | 8 | data.frame | list |
| weights | 14 | -none- | numeric |
| k2prime | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(mrtm2out$fit), digits = 3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 4.86 | 27622 | 4.41e-21 | 2 | -39.2 | 84.4 |
| BIC | deviance | df.residual |
|---|---|---|
| 86.1 | 260 | 11 |
\[\int _{0} ^{t} C_{T} (\tau) d\tau = DVR\Bigg(\int _{0} ^{t} C_{T}' (\tau) d\tau+\frac{C_{T}'(t)}{k2'}\Bigg)-bC_{T}(t)\]
refmlLogan is the kinetic model used in WAPI, which is detailed in Turkheimer’s 2003 IEEE paper. This model requires input of the times, reftac and roitac and k2prime. Weights are optional.
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
taclow = tacdata$FSLSFC # Frontal Cortex
tacmed = tacdata$fslLST # Limbic Striatum
tachigh = tacdata$fslAST # Associative Striatum
refmlLogan_tstar(times, reftac, taclow, tacmed, tachigh, k2prime, 'demonstration')Using this information (unshown), we see that we should probably use 9 included frames for the fit. We shall use 6.
refmlloganout <- refmlLogan(times, reftac, roitac, k2prime,
tstarIncludedFrames = 9, weights=weights)pandoc.table(refmlloganout$par, digits = 3)| bp | k2 |
|---|---|
| 1.53 | 0.0437 |
The output object can be used to visualise the fit:
plot_refmlLoganfit(refmlloganout)RefmlLogan Fit Check
… and the residuals:
plot_residuals(refmlloganout)RefLogan Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(refmlloganout))| Length | Class | Mode | |
|---|---|---|---|
| par | 2 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 3 | data.frame | list |
| fitvals | 4 | data.frame | list |
| weights | 14 | -none- | numeric |
| k2prime | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(refmlloganout$fit), digits = 3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 94.4 | 48778 | 3.13e-15 | 2 | -53.2 | 112 |
| BIC | deviance | df.residual |
|---|---|---|
| 113 | 62431 | 7 |
refPatlak requires input of the times, reftac and roitac. Weights are optional.
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
taclow = tacdata$FSLSFC # Frontal Cortex
tacmed = tacdata$fslLST # Limbic Striatum
tachigh = tacdata$fslAST # Associative Striatum
refPatlak_tstar(times, reftac, taclow, tacmed, tachigh, 'demonstration')Since this is actually a reversible tracer, we shouldn’t really be using refPatlak at all, but we can use four frames to test it out (unshown).
refpatlakout <- refPatlak(times, reftac, roitac,
tstarIncludedFrames = 4, weights=weights)pandoc.table(refpatlakout$par, digits=4)##
## --------
## K
## --------
## 0.005227
## --------
The output object can be used to visualise the fit:
plot_refPatlakfit(refpatlakout)RefPatlak Fit Check
… and the residuals:
plot_residuals(refpatlakout)RefPatlak Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(refpatlakout))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 3 | data.frame | list |
| fitvals | 2 | data.frame | list |
| weights | 14 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(refpatlakout$fit), digits = 3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 0.894 | 0.84 | 0.0435 | 16.8 | 0.0547 | 2 | 8.06 | -10.1 |
| BIC | deviance | df.residual |
|---|---|---|
| -12 | 0.00378 | 2 |
First all the data is loaded (unshown). This data includes:
The data is plotted below:
par(mfrow=c(2,2))
plot(t_tac, tac)
plot(blooddata$Time.sec., blooddata$Cbl.nCi.cc., xlab = 'Time (sec)', ylab = 'Cb')
plot(plasmadata$Time.sec., plasmadata$Cpl.nCi.cc., xlab = 'Time (sec)', ylab = 'Cp')
plot(parentdata$Times, parentdata$Fraction, xlab = 'Time (sec)', ylab = 'Parent Fraction')Compartmental Model Data
First, we must organise the blood data and put it together. This function can also accept parent-fraction-corrected plasma data too.
input <- blood_interp(t_blood = blooddata$Time.sec./60,
blood = blooddata$Cbl.nCi.cc.,
t_plasma = plasmadata$Time.sec./60,
plasma = plasmadata$Cpl.nCi.cc.,
t_parentfrac = parentdata$Times/60,
parentfrac = parentdata$Fraction)1TCM requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function and vB are also optional. If they are not included, then the function will fit them.
onetcmout <- onetcm(t_tac, tac, input, weights=weights)pandoc.table(onetcmout$par, digits = 3)| K1 | k2 | inpshift | vB | Vt |
|---|---|---|---|---|
| 0.112 | 0.0317 | -0.164 | 0.0824 | 3.54 |
The inpshift is the fitted delay. We can visualise the matching of the two curves without any time shifting (i.e. as they are in the same time scale):
plot_inptac_timings(t_tac, tac, input, inpshift = 0)Delay Check Before
The timing of these two clearly do not match. Now we can visualise the fit with the time shifting fitted by the one tissue compartment model: in this case, the delay is negative, which means that the TAC was shifted to be a little bit later.
plot_inptac_timings(t_tac, tac, input, inpshift = onetcmout$par$inpshift)Delay Check After
This parameter has clearly performed a good job in fitting the timing of the two curves to one another.
The output object can be used to visualise the model fit:
plot_1tcmfit(onetcmout)1TCM Fit Check
… and the residuals:
plot_residuals(onetcmout)1TCM Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(onetcmout))| Length | Class | Mode | |
|---|---|---|---|
| par | 5 | data.frame | list |
| par.se | 5 | data.frame | list |
| fit | 6 | nls | list |
| tacs | 3 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift_fitted | 1 | -none- | logical |
| vB_fitted | 1 | -none- | logical |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(onetcmout$fit), digits=3)| sigma | isConv | finTol | logLik | AIC | BIC | deviance | df.residual |
|---|---|---|---|---|---|---|---|
| 22.2 | TRUE | 1.49e-08 | -182 | 373 | 381 | 14759 | 30 |
2TCM requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function and vB are also optional. If they are not included, then the function will fit them. In this case, we can use the delay from the 1tcm for the 2tcm fit, and it will not be fitted.
twotcmout <- twotcm(t_tac, tac, input, weights=weights, inpshift = onetcmout$par$inpshift)pandoc.table(twotcmout$par, digits = 3)| K1 | k2 | k3 | k4 | vB | Vt |
|---|---|---|---|---|---|
| 0.128 | 0.0721 | 0.0436 | 0.0295 | 0.0503 | 4.42 |
The output object can be used to visualise the model fit:
plot_2tcmfit(twotcmout)2TCM Fit Check
… and the residuals:
plot_residuals(twotcmout)2TCM Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(twotcmout))| Length | Class | Mode | |
|---|---|---|---|
| par | 6 | data.frame | list |
| par.se | 6 | data.frame | list |
| fit | 6 | nls | list |
| tacs | 3 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift_fitted | 1 | -none- | logical |
| vB_fitted | 1 | -none- | logical |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(twotcmout$fit), digits=3)| sigma | isConv | finTol | logLik | AIC | BIC | deviance | df.residual |
|---|---|---|---|---|---|---|---|
| 7.73 | TRUE | 1.49e-08 | -141 | 294 | 303 | 1731 | 29 |
The Logan Plot method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function and vB are also optional. If they are not included, then they will be set to zero. For the input delay, we can use the delay from the 1tcm fit. If vB is included, it will be corrected for before the parameter estimation using the following formula:
\[ C_{T}(t) = \frac{C_{Measured}(t) - vB\times C_{B}(t)}{1-vB} \]
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
lowroi = inptacdata$PUT
medroi = inptacdata$GM
highroi = inptacdata$gmCER
Logan_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration',
inpshift = onetcmout$par$inpshift, vB = 0.05, gridbreaks=4)Using this information (unshown), we see that we should probably use 10 frames for the fit.
loganout <- Loganplot(t_tac, tac, input, 10, weights, inpshift = onetcmout$par$inpshift, vB=0.05)loganout$par## Vt
## 1 4.225652
The output object can be used to visualise the model fit:
plot_Loganfit(loganout)Logan Plot Fit Check
… and the residuals:
plot_residuals(loganout)Logan Plot Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(loganout))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| par.se | 1 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 2 | data.frame | list |
| fitvals | 2 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift | 1 | -none- | numeric |
| vB | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(loganout$fit), digits=3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 0.998 | 0.998 | 1.49 | 4754 | 2.18e-12 | 2 | -17.2 | 40.3 |
| BIC | deviance | df.residual |
|---|---|---|
| 41.2 | 17.8 | 8 |
The MA1 method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional. If it is not included, then it will be set to zero.
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
lowroi = inptacdata$PUT
medroi = inptacdata$GM
highroi = inptacdata$gmCER
ma1_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration',
inpshift = onetcmout$par$inpshift, gridbreaks=4)## png
## 2
Using this information, we see that we should probably use very few frames for the fit: 11 seems reasonably good and reasonably stable.
ma1out <- ma1(t_tac, tac, input, 11, weights, inpshift = onetcmout$par$inpshift)ma1out$par## Vt
## 1 4.143303
The output object can be used to visualise the model fit:
plot_ma1fit(ma1out)MA1 Plot Fit Check
… and the residuals:
plot_residuals(ma1out)MA1 Plot Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(ma1out))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 2 | data.frame | list |
| fitvals | 6 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(ma1out$fit), digits=3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 0.999 | 0.999 | 6.6 | 4202 | 4.28e-14 | 2 | -35.4 | 76.7 |
| BIC | deviance | df.residual |
|---|---|---|
| 77.9 | 392 | 9 |
The MA2 method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional. If it is not included, then it will be set to zero.
ma2out <- ma2(t_tac, tac, input, weights, inpshift = onetcmout$par$inpshift)pandoc.table(ma2out$par, digits = 3)| Vt | Vs | Vnd | K1 | k2 | k3 | k4 |
|---|---|---|---|---|---|---|
| 4.11 | 3.42 | 0.684 | 0.153 | 0.223 | 0.00187 | 0.000373 |
One can also calculate the rate constants from the fit terms. The article provides a formula for calculation of VS, but they left out a number though for which term to use (this error is also on the PMOD page, which simply displays an image of the formula from the article). I believe I have chosen the correct one, but this might not be the case. VND can be calculated by subtracting Vs from VT.
The output object can be used to visualise the model fit:
plot_ma2fit(ma2out)MA2 Plot Fit Check
… and the residuals:
plot_residuals(ma2out)MA2 Plot Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(ma2out))| Length | Class | Mode | |
|---|---|---|---|
| par | 7 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 2 | data.frame | list |
| fitvals | 8 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(ma2out$fit), digits=3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 0.977 | 0.974 | 40.6 | 319 | 4.09e-24 | 4 | -172 | 355 |
| BIC | deviance | df.residual |
|---|---|---|
| 363 | 49558 | 30 |
\[\int _{0} ^{t} C_{T} (\tau) d\tau = V_{T}\int _{0} ^{t} C_{P} (\tau) d\tau-\frac{1}{K_{2}}C_{T}(t)\]
The multilinear Logan Plot is the kinetic model used in WAPI, which is detailed in Turkheimer’s 2003 IEEE paper. This model requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional, and if it is not included, it will be set to zero.
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
lowroi = inptacdata$PUT
medroi = inptacdata$GM
highroi = inptacdata$gmCER
mlLogan_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration',
inpshift = onetcmout$par$inpshift, gridbreaks=4)Using this information (unshown), we see that we should probably use 10 frames for the fit.
mlloganout <- mlLoganplot(t_tac, tac, input, 10, weights, inpshift = onetcmout$par$inpshift)pandoc.table(mlloganout$par, digits=3)| Vt | k2 |
|---|---|
| 4.21 | 0.0179 |
The output object can be used to visualise the model fit:
plot_mlLoganfit(mlloganout)mlLogan Plot Fit Check
… and the residuals:
plot_residuals(mlloganout)mlLogan Plot Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(mlloganout))| Length | Class | Mode | |
|---|---|---|---|
| par | 2 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 2 | data.frame | list |
| fitvals | 4 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
| model | 1 | -none- | character |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(mlloganout$fit), digits=3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 278 | 14762 | 5.39e-15 | 2 | -69.4 | 145 |
| BIC | deviance | df.residual |
|---|---|---|
| 146 | 617196 | 8 |
The Patlak Plot method requires input of the TAC times, TAC and input. Weights are optional. The time shift of the input function is also optional. If it is not included, then it will be set to zero. For the input delay, we can use the delay from the 1tcm fit.
First, we wish to ascertain what a suitable t* value would be. For this, we use TACs from a low-, medium- and high-binding region.
Patlak_tstar(t_tac, lowroi, medroi, highroi, input, filename='demonstration',
inpshift = onetcmout$par$inpshift, gridbreaks=4)Since this is actually a reversible tracer, we shouldn’t really be using the Patlak Plot at all, but we can use nine frames to test it out (unshown).
patlakout <- Patlakplot(t_tac, tac, input, 9, weights, inpshift = onetcmout$par$inpshift)pandoc.table(patlakout$par, digits=3)| K |
|---|
| 0.0118 |
The output object can be used to visualise the model fit:
plot_Patlakfit(patlakout)Patlak Plot Fit Check
… and the residuals:
plot_residuals(patlakout)Patlak Plot Residuals
Included in the output are all the parameters, the parameters of the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(patlakout))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| par.se | 1 | data.frame | list |
| fit | 13 | lm | list |
| tacs | 2 | data.frame | list |
| fitvals | 2 | data.frame | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| inpshift | 1 | -none- | numeric |
| tstarIncludedFrames | 1 | -none- | numeric |
This means that the fit can be easily evaluated and compared with other fits:
pandoc.table(glance(patlakout$fit), digits=3)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC |
|---|---|---|---|---|---|---|---|
| 0.997 | 0.997 | 0.0828 | 2772 | 2.34e-10 | 2 | 10.7 | -15.4 |
| BIC | deviance | df.residual |
|---|---|---|
| -14.8 | 0.048 | 7 |
The SIME method requires input of the TAC times, several TACs, input and a grid of VND values. Weights are optional. ROI sizes can optionally be included to weight the contributions of the different ROIs. The time shift of the input function and vB are also optional. If they are not included, they will be set to 0 and 0.05 respectively.
tacdf <- data.frame(GM = inptacdata$GM,
ACC = inptacdata$gmACC,
CER = inptacdata$gmCER,
PUT = inptacdata$PUT)
Vndgrid <- seq(from=0, to=5, by=0.1)
SIMEout <- SIME(t_tac, tacdf, input, Vndgrid, weights = weights,
inpshift = onetcmout$par$inpshift, vB = 0.05, twotcmstart = 1)pandoc.table(SIMEout$par, digits=3)| Vnd |
|---|
| 1.8 |
The output object can be used to visualise the model fit:
plot_SIMEfit(SIMEout)SIME Plot Fit Check
Included in the output are all the parameters, the parameters used for the fitting functions, the fits themselves and the raw data.
pandoc.table(summary(SIMEout))| Length | Class | Mode | |
|---|---|---|---|
| par | 1 | data.frame | list |
| tacs | 7 | data.frame | list |
| fitvals | 2 | data.frame | list |
| roifits | 5 | tbl_df | list |
| input | 4 | data.frame | list |
| weights | 38 | -none- | numeric |
| roiweights | 4 | -none- | numeric |
| inpshift | 1 | -none- | numeric |
| vB | 1 | -none- | numeric |
| model | 1 | -none- | character |
Alternatively, this model can be implemented in a faster manner of sequentially refining its results:
Vndgrid <- seq(from=0, to=5, by=1)
SIMEout1 <- SIME(t_tac, tacdf, input, Vndgrid, weights = weights,
inpshift = onetcmout$par$inpshift)
Vndgrid <- seq(from=SIMEout1$par$Vnd-1, to=SIMEout1$par$Vnd+1, length.out = 10)
SIMEout2 <- SIME(t_tac, tacdf, input, Vndgrid, weights = weights,
inpshift = onetcmout$par$inpshift)
Vndgrid <- seq(from=SIMEout2$par$Vnd-0.2, to=SIMEout2$par$Vnd+0.2, length.out = 10)
SIMEout3 <- SIME(t_tac, tacdf, input, Vndgrid, weights = weights,
inpshift = onetcmout$par$inpshift)